#include "classes.h"

#define EM_E0 GSL_CONST_MKSA_VACUUM_PERMITTIVITY //8.854187817e-12
#define EM_M0 GSL_CONST_MKSA_VACUUM_PERMEABILITY //1.25663706144e-6

double Integrals::cosa(double fis)
{
	double val;
	val=sin(fi-fis)*sin(teta);
	return val;
}

double Integrals::cosb(double fis)
{
	double val;
	val=sin(fi-fis)*cos(teta);
	return val;
}

double Integrals::cosg(double fis)
{
	double val;
	val=cos(fi-fis)*sin(teta);
	return val;
}

double Integrals::cosd(double fis)
{
	double val;
	val=cos(fi-fis)*cos(teta);
	return val;
}

__complex__ double Integrals::ji1(double ro)
{
	__complex__ double val;
	val=1I*gsl_sf_bessel_I1(__imag__(alfa2)*ro);
	return val;
}

__complex__ double Integrals::yi1(double ro)
{
	__complex__ double val;
        val=-1*gsl_sf_bessel_I1(__imag__(alfa2)*ro)+2I/M_PI*gsl_sf_bessel_K1(__imag__(alfa2)*ro);
	return val;
}

double Integrals::js1(double ro)
{
	double val;
	val=gsl_sf_bessel_J0(alfa1*ro)-gsl_sf_bessel_J1(alfa1*ro)/alfa1/ro;
	return val;
}

double Integrals::ys1(double ro)
{
	double val;
	val=gsl_sf_bessel_Y0(alfa1*ro)-gsl_sf_bessel_Y1(alfa1*ro)/alfa1/ro;
	return val;
}

double Integrals::js2(double ro)
{
	double val;
	val=gsl_sf_bessel_I0(__imag__(alfa2)*ro)-gsl_sf_bessel_I1(__imag__(alfa2)*ro)/__imag__(alfa2)/ro;
	return val;
}

__complex__ double Integrals::ys2(double ro)
{
	__complex__ double val,y0;
        y0=1I*gsl_sf_bessel_I0(__imag__(alfa2)*ro)-2/M_PI*gsl_sf_bessel_K0(__imag__(alfa2)*ro);
	val=y0-yi1(ro)/alfa2/ro;
	return val;
}

__complex__ double Integrals::cexp(double rr)
{
	__complex__ double val=1I;
	val=cos(k*rr)-val*sin(k*rr);
	val/=rr;
	return val;
}

__complex__ double Integrals::cexp0(double ro, double fis)
{
	__complex__ double val=1I;
	val=cos(k*ro*cosg(fis))+val*sin(k*ro*cosg(fis));
	return val;
}

__complex__ double Integrals::G(int n)
{
	__complex__ double val;
	float ee;
	if(n==1) ee=e1; else ee=e2;
	val=ee*js1(a)/gsl_sf_bessel_J1(alfa1*a)/alfa1-Rs1(a)/R1(a)/alfa2;
        val*=w*EM_E0*a/beta/(1/pow(alfa1,2)+1/pow(__imag__(alfa2),2));	//Только для чистомнимых alfa2
	return val;
}

double Integrals::R1(double ro)
{
	double val;
        I1r=gsl_sf_bessel_I1(__imag__(alfa2)*ro);
	K1r=gsl_sf_bessel_K1(__imag__(alfa2)*ro);
        val=2/M_PI*(I1b*K1r-I1r*K1b);
	return val;
}

__complex__ double Integrals::Rs1(double ro)
{
	__complex__ double val;
        I1r=gsl_sf_bessel_I1(__imag__(alfa2)*ro);
	K1r=gsl_sf_bessel_K1(__imag__(alfa2)*ro);
	I0r=gsl_sf_bessel_I0(__imag__(alfa2)*ro);
	K0r=gsl_sf_bessel_K0(__imag__(alfa2)*ro);
	ro=__imag__(alfa2)*ro;
        val=2I/M_PI*(I0r*K1b-I1r*K1b/ro+I1b*K0r+I1b*K1r/ro);
	return val;
}

__complex__ double Integrals::S1(double ro)
{
	__complex__ double val;
        I1r=gsl_sf_bessel_I1(__imag__(alfa2)*ro);
	K1r=gsl_sf_bessel_K1(__imag__(alfa2)*ro);
	I0r=gsl_sf_bessel_I0(__imag__(alfa2)*ro);
	K0r=gsl_sf_bessel_K0(__imag__(alfa2)*ro);
	ro=__imag__(alfa2)*b;
        val=2I/M_PI*(I1r*K1b/ro-I1r*K0b+I1b*K1r/ro-I0b*K1r);
	return val;
}

double Integrals::Ss1(double ro)
{
	double val;
	double rob;
	rob=__imag__(alfa2)*b;
        I1r=gsl_sf_bessel_I1(__imag__(alfa2)*ro);
	K1r=gsl_sf_bessel_K1(__imag__(alfa2)*ro);
	I0r=gsl_sf_bessel_I0(__imag__(alfa2)*ro);
	K0r=gsl_sf_bessel_K0(__imag__(alfa2)*ro);
	ro=__imag__(alfa2)*ro;
        val=2/M_PI*((I0b-I1b/rob)*(K0r+K1r/ro)-(I0r-I1r/ro)*(K0b+K1b/rob));
	return val;
}

__complex__ double Integrals::Ero1(double ro,double fis)
{
	__complex__ double val;
        val=1I*(-beta*alfa1*js1(ro)+w*EM_M0*G1*gsl_sf_bessel_J1(alfa1*ro)/ro)*sin(fis);
	return val;	
}

__complex__ double Integrals::Efi1(double ro, double fis)
{
	__complex__ double val;
        val=1I*(-beta*gsl_sf_bessel_J1(alfa1*ro)/ro+w*EM_M0*G1*alfa1*js1(ro))*cos(fis);
	return val;	
}

__complex__ double Integrals::Hro1(double ro, double fis)
{
	__complex__ double val;
        val=1I*(-beta*G1*alfa1*js1(ro)+w*EM_E0*e1*gsl_sf_bessel_J1(alfa1*ro)/ro)*cos(fis);
	return val;	
}

__complex__ double Integrals::Hfi1(double ro, double fis)
{
	__complex__ double val;
        val=1I*(beta*G1*gsl_sf_bessel_J1(alfa1*ro)/ro-w*EM_E0*e1*alfa1*js1(ro))*sin(fis);
	return val;	
}

__complex__ double Integrals::Ero2(double ro,double fis)
{
	__complex__ double val;
        val=1I*pow(alfa1,2)/alfa2*gsl_sf_bessel_J1(alfa1*a)*(-beta/R1(a)*Rs1(ro)+w*EM_M0*G2/ro/alfa2/S1(a)*S1(ro))*sin(fis);
	return val;	
}

__complex__ double Integrals::Efi2(double ro, double fis)
{
	__complex__ double val;
        val=1I*pow(alfa1,2)/alfa2*gsl_sf_bessel_J1(alfa1*a)*(-beta/R1(a)/alfa2*R1(ro)/ro+w*EM_M0*G2/S1(a)*Ss1(ro))*cos(fis);
	return val;	
}

__complex__ double Integrals::Hro2(double ro, double fis)
{
	__complex__ double val;
        val=1I*pow(alfa1,2)/alfa2*gsl_sf_bessel_J1(alfa1*a)*(-beta*G2/S1(a)*Ss1(ro)+w*EM_E0*e2/alfa2/R1(a)*R1(ro)/ro)*cos(fis);
	return val;	
}

__complex__ double Integrals::Hfi2(double ro, double fis)
{
	__complex__ double val;
        val=1I*pow(alfa1,2)/alfa2*gsl_sf_bessel_J1(alfa1*a)*(beta*G2/alfa2/S1(a)*S1(ro)/ro-w*EM_E0*e2/R1(a)*Rs1(ro))*sin(fis);
	return val;	
}

double Integrals::alfa()
{
	int j=0,q;
	double alf=1000,xa,xb,det,res=0;
        __complex__ double d1,d2;
        det=10;
	while (j<1)
	{
		xa=alf;	xb=xa+det;
		//first point
		alfa1=xa;
		beta=sqrt(pow(k,2)*e1*m1-pow(alfa1,2));
                alfa2=sqrt(pow(beta,2)-pow(k,2)*e2*m2)+0I;
		alfa2=1I*alfa2;
		I1b=gsl_sf_bessel_I1(__imag__(alfa2)*b);
		I0b=gsl_sf_bessel_I0(__imag__(alfa2)*b);
		K1b=gsl_sf_bessel_K1(__imag__(alfa2)*b);
		K0b=gsl_sf_bessel_K0(__imag__(alfa2)*b);
		d1=(alfa1*a*js1(a)/gsl_sf_bessel_J1(alfa1*a)-m2/m1*pow(alfa1*a,2)/alfa2/a*Ss1(a)/S1(a));
		d1=d1*(alfa1*a*js1(a)/gsl_sf_bessel_J1(alfa1*a)-e2/e1*pow(alfa1*a,2)/alfa2/a*Rs1(a)/R1(a));
		d1=d1-pow(beta/k,2)*pow(1+pow(alfa1,2)/pow(__imag__(alfa2),2),2)/e1/m1;
		//second point
		alfa1=xb;
		beta=sqrt(pow(k,2)*e1*m1-pow(alfa1,2));
                alfa2=sqrt(pow(beta,2)-pow(k,2)*e2*m2)+0I;
		alfa2=1I*alfa2;
		I1b=gsl_sf_bessel_I1(__imag__(alfa2)*b);
		I0b=gsl_sf_bessel_I0(__imag__(alfa2)*b);
		K1b=gsl_sf_bessel_K1(__imag__(alfa2)*b);
		K0b=gsl_sf_bessel_K0(__imag__(alfa2)*b);
		d2=alfa1*a*js1(a)/gsl_sf_bessel_J1(alfa1*a)-m2/m1*pow(alfa1*a,2)/alfa2/a*Ss1(a)/S1(a);
		d2=d2*(alfa1*a*js1(a)/gsl_sf_bessel_J1(alfa1*a)-e2/e1*pow(alfa1*a,2)/alfa2/a*Rs1(a)/R1(a));
		d2=d2-pow(beta/k,2)*pow(1+pow(alfa1,2)/pow(__imag__(alfa2),2),2)/e1/m1;
		if (GSL_SIGN(__real__(d1))!=GSL_SIGN(__real__(d2)))
		{
			for (q=0;q<16;q++)
			{
				alf=(xa+xb)/2;
				// first point
				alfa1=xa;
				beta=sqrt(pow(k,2)*e1*m1-pow(alfa1,2));
                                alfa2=sqrt(pow(beta,2)-pow(k,2)*e2*m2)+0I;
				alfa2=1I*alfa2;
				I1b=gsl_sf_bessel_I1(__imag__(alfa2)*b);
				I0b=gsl_sf_bessel_I0(__imag__(alfa2)*b);
				K1b=gsl_sf_bessel_K1(__imag__(alfa2)*b);
				K0b=gsl_sf_bessel_K0(__imag__(alfa2)*b);
				d1=(alfa1*a*js1(a)/gsl_sf_bessel_J1(alfa1*a)-m2/m1*pow(alfa1*a,2)/alfa2/a*Ss1(a)/S1(a));
				d1=d1*(alfa1*a*js1(a)/gsl_sf_bessel_J1(alfa1*a)-e2/e1*pow(alfa1*a,2)/alfa2/a*Rs1(a)/R1(a));
				d1=d1-pow(beta/k,2)*pow(1+pow(alfa1,2)/pow(__imag__(alfa2),2),2)/e1/m1;
				//second point
				alfa1=alf;
				beta=sqrt(pow(k,2)*e1*m1-pow(alfa1,2));
                                alfa2=sqrt(pow(beta,2)-pow(k,2)*e2*m2)+0I;
				alfa2=1I*alfa2;
				I1b=gsl_sf_bessel_I1(__imag__(alfa2)*b);
				I0b=gsl_sf_bessel_I0(__imag__(alfa2)*b);
				K1b=gsl_sf_bessel_K1(__imag__(alfa2)*b);
				K0b=gsl_sf_bessel_K0(__imag__(alfa2)*b);
				d2=(alfa1*a*js1(a)/gsl_sf_bessel_J1(alfa1*a)-m2/m1*pow(alfa1*a,2)/alfa2/a*Ss1(a)/S1(a));
				d2=d2*(alfa1*a*js1(a)/gsl_sf_bessel_J1(alfa1*a)-e2/e1*pow(alfa1*a,2)/alfa2/a*Rs1(a)/R1(a));
				d2=d2-pow(beta/k,2)*pow(1+pow(alfa1,2)/pow(__imag__(alfa2),2),2)/e1/m1;
				if (GSL_SIGN(__real__(d1))!=GSL_SIGN(__real__(d2)))	xb=alf;
				else	xa=alf;
			}
			res=alf;
			j++;
                }
		alf=alf+det;
	}
        return	res;
}
